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The development of hyperbolic formulations of Einstein’s equations has revolutionized our ability 
to perform long-time, stable, accurate numerical simulations of strong field gravitational phenomena. 
However, hyperbolic methods have seen relatively little application in one area of interest, type II 
critical collapse, where the challenges for a numerical code are particularly severe. Using the critical 
collapse of a massless scalar field in spherical symmetry as a test case, we study a generalization of 
the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formulation due to Brown that is suited for use 
with curvilinear coordinates. We adopt standard dynamical gauge choices, including H-log slicing 
and a shift that is either zero or evolved by a Gamma-driver condition. With both choices of shift 
we are able to evolve sufficiently close to the black hole threshold to (1) unambiguously identify 
the discrete self-similarity of the critical solution, (2) determine an echoing exponent consistent 
with previous calculations, and (3) measure a mass scaling exponent, also in accord with prior 
computations. Our results can be viewed as an encouraging first step towards the use of hyperbolic 
formulations in more generic type II scenarios, including the as yet unresolved problem of critical 
collapse of axisymmetric gravitational waves. 

PACS numbers: 04.25.dc, 04.40.-b, 04.40.Dg 


I. INTRODUCTION 

In this paper we investigate the application of the 
Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formula¬ 
tion of Einstein’s equations [10, as well as the dynami¬ 
cal coordinate choices typically associated with it, within 
the context of critical gravitational collapse. The BSSN 
formulation is a recastiM of the standard 3-1-1 Arnowitt- 
Deser-Misner (ADlVn H equations that is known to be 
strongly hyperbolic [J, and suitable for numerical stud¬ 
ies. It has been widely used in numerical relativity and 
provides a robust and stable evolution for the spacetime 
geometry. Most notably, various implementations of this 
formulation have allowed successful computation of dy¬ 
namical spacetimes describing binaries of gravitationally- 
compact objects @-0. The standard gauge choices in 
BSSN—namely the l-blog slicing condition 0 and the 
Gamma-driver shift condition [l0|— are partial differen¬ 
tial equations (PDEs) of evolutionary type. Furthermore, 
the BSSN approach results in a set of so-called free evo¬ 
lution equations, meaning that the Hamiltonian and mo¬ 
mentum constraints are only solved at the initial time. 
Thus, once initial data has been determined, one only 
has to solve time-dependent PDEs in order to compute 
the geometric variables in the BSSN scheme. In partic¬ 
ular, during the evolution there is no need to solve any 
elliptic equations, which in general could arise either from 
the constraints or from coordinate conditions. This is ad¬ 
vantageous since it can be quite challenging to implement 
efficient numerical elliptic solvers. 

In addition to the BSSN approach, the numerical rela¬ 
tivity community has adopted the generalized harmonic 
(GH) [HI formulation of Einstein’s equations, which is 
also strongly hyperbolic and has performed very well in 
simulations of compact binaries [1, [i2|- Like BSSN, the 


GH formulation is of evolutionary type so that all of the 
metric components satisfy time-dependent PDEs. It too 
uses dynamical coordinate choices: in this case one needs 
to provide a prescription for the evolution of the har¬ 
monic functions defined by 

Despite the tremendous success of these hyperbolic 
formulations in evolving strongly gravitating spacetimes 
containing black holes and neutron stars, they have not 
seen widespread use in another area of strong gravity 
physics typically studied via numerical relativity, namely 
critical phenomena in gravitational collapse. First re¬ 
ported in and briefly reviewed below, critical phe¬ 
nomena emerge at the threshold of black hole formation 
and present significant challenges for thorough and accu¬ 
rate computational treatment. The original observation 
of critical behaviour as well as many of the subsequent 
studies were restricted to spherical symmetry (for a re¬ 
view, see [13, [3) and there is a clear need to extend the 
work to more generic cases. In this respect the BSSN and 
GH formulations would appear to be attractive frame¬ 
works. However, it is not yet clear if these hyperbolic 
formulations, in conjunction with the standard dynami¬ 
cal gauge choices that have been developed, will allow the 
critical regime to be probed without the development of 
coordinate pathologies. Particularly notable in this re¬ 
gard is an implementation of the GH formulation that 
was employed by Sorkin and Choptuik [l^ to study the 
critical collapse of a massless scalar field in spherical sym¬ 
metry. Despite extensive experimentation with a variety 
of coordinate conditions, the code that was developed was 
not able to calculate near-critical spacetimes: coordinate 
singularities invariably formed once the critical regime 
was approached. A natural question that then arises is 
whether the BSSN formulation (including the standard 
dynamical gauge choices used with it) is similarly prob- 
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lematic or if it provides an effective framework to study 
critical phenomena. 

Here we begin the task of addressing this question by 
revisiting the model of spherically symmetric massless 
scalar collapse. We use a generalization of the BSSN 
formulation due to Brown [l^ that is well suited for use 
with curvilinear coordinates. The choice of a massless 
scalar field as the matter source has the great advantage 
that the nature of the critical solution is very well known 
[nil, making it straightforward for us to determine if 
and when our approach has been successful. We note that 
although the calculations described below are restricted 
to spherical symmetry our ultimate goal is to develop an 
evolutionary scheme—including gauge choices—that can 
be applied to a variety of critical phenomena studies in 
axial symmetry and ultimately generic cases. 

We now briefly review the main concepts and features 
of black hole critical phenomena that are most pertinent 
to the work in this paper. Full details and pointers to 
the extensive literature on the subject may be found in 
review articles [l3, [III . 

Critical phenomena in gravitational collapse can be 
described as a phase transition, analogous to that in a 
thermodynamical system. Under certain assumptions, a 
matter source coupled to the Einstein gravitational field 
will evolve to one of two distinct final phases. On the one 
hand, weak initial data will eventually disperse to infin¬ 
ity leaving flat spacetime as the end state. On the other 
hand, sufficiently strong data will develop significant self 
gravitation and then collapse, resulting in a final phase 
which contains a black hole. Quite generically, remark¬ 
able behaviour emerges at and near the transition be¬ 
tween these phases, and this behaviour is precisely what 
we mean by the critical phenomena in the system under 
consideration. 

It transpires that there are two broad classes of critical 
phenomena that can be distinguished by the behaviour 
of the black hole mass at threshold. The class of interest 
here, known as type II, is characterized by infinitesimal 
mass at the transition. Further, the black hole mass, 
Mbh, satisfies a scaling law: 

Mbh'^\p-P*P, ( 1 ) 

where p is an arbitrary parameter that controls the 
strength of the matter source at the initial time, p* is the 
parameter value at threshold and the mass scaling expo¬ 
nent, 7, is a constant that is independent of the choice of 
the initial data. Type II behaviour is also characterized 
by the emergence of a unique solution at threshold which 
is generically self-similar. In some cases, including the 
massless scalar field, the self-similarity is discrete. Specif¬ 
ically, in spherically symmetric critical collapse with dis¬ 
crete self-similarity (DSS), as p —>■ p* we find 

Z*{p + A,t + A)^Z*{p,t), (2) 

where Z* represents some scale-invariant component 
(function) of the critical solution. Here p = In(rs) 


and T = ln(rs ~ ^s) logarithmically rescaled values 
of the areal radius, rs, and polar time, Tg, respectively, 
and Tg is the accumulation time at which the central sin¬ 
gularity associated with the DSS solution forms. Tg has 
been normalized so that it measures proper time at the 
origin. As with 7, the echoing (rescaling) exponent, A, 
is a universal constant for a specific matter source; i.e. it 
is independent of the form of the initial data. 

Another feature of type H collapse, intimately related 
to the self-similarity of the critical solution, is that the 
curvature can become arbitrarily large: in the limit of 
infinite fine-tuning, p —>■ p*, a naked singularity forms. 
Furthermore, the echoing behaviour @ results in the 
development of fine structure in the solution around the 
center of the scaling symmetry. Observing this structure 
and measuring the echoing exponent A associated with 
it requires a code that can reliably evolve solutions very 
close to the critical spacetime and that provides sufficient 
numerical resolution in the vicinity of the accumulation 
point (rs,Ts) = (0,rg*). 

As mentioned above, most studies of critical phenom¬ 
ena have assumed spherical symmetry. This is particu¬ 
larly so for the case of type H behaviour where the resolu¬ 
tion demands dictated by the self-similarity of the critical 
solutions makes multi-dimensional work extremely com¬ 
putationally intensive. As far as we know, the only work 
in spherical symmetry to have used a purely evolutionary 
approach based on the BSSN or GH forms of the Einstein 
equations is [l^ which, as we have noted, was not suc¬ 
cessful in isolating the critical solution.^ In axisymme- 
try there have been two investigations of type H collapse 
of massless scalar fields [2^, [2^ , and several of type H 
collapse of pure gravitational waves (vacuum) [27H3lll . 
Of these, only Alcubierre et al.’s [2^ and Sorkin’s |3l| 
calculations of vacuum collapse adopted hyperbolic for¬ 
malisms, and only the scalar field calculations—which 
employed a modified ADM formulation and partially con¬ 
strained evolution—were able to completely resolve the 
critical behaviour, including the discrete self-similarity of 
the critical solutions. In the fully three-space dimension 
(3D) context there have also been a few studies of type 
H collapse to date. Perhaps most notable is the recent 
work of Healy and Laguna which used a massless 
scalar field as a matter source and the BSSN formulation 
with standard dynamical gauge choices. The authors 
were able to observe the mass scaling © with a mea¬ 
sured 7 « 0.37 consistent with calculations in spherical 
symmetry. However, they were not able to conclusively 
see the discrete self-similarity of the critical solution; in 
particular they could not accurately measure the echo¬ 
ing exponent, A. This shortcoming was attributed to a 
lack of computational resources rather than a breakdown 


^ However, see l24l for an investigation of type II behaviour in the 
collapse of a scalar field in 2-1-1 AdS spacetime that employs an 
ad hoc free evolution scheme. 



3 


of the underlying methodology, including the coordinate 
conditions that were adopted. Finally, there have been 
two attempts to probe the black hole threshold for the 
collapse of pure gravitational waves in 3D [s^ . Both 

employed a BSSN approach with, for the most part, stan¬ 
dard dynamical gauge choices. In both cases problems 
with the gauge apparently precluded calculation near the 
critical point (although resolution limitations may also 
have been an issue) and neither the mass scaling nor the 
echoing exponent could be be estimated in either study. 

We can thus summarize the state of the art in the use 
of hyperbolic formulations for the study of type II criti¬ 
cal collapse as follows: to our knowledge there has been 
no implementation of a fully evolutionary scheme, based 
on either BSSN or GH, that has allowed for evolution 
sufficiently close to a precisely critical solution to allow 
the unambiguous identification of discrete self-similarity 
(or continuous self-symmetry for that matter). Again, 
and particularly in light of the experience of [l6| , the key 
aim of this paper is to investigate the extent to which it 
is possible to use a BSSN scheme to fully resolve type II 
solutions. A major concern here is the appropriate choice 
of coordinate conditions, not least since dynamical gauge 
choices can be prone to the development of gauge shocks 
and other types of coordinate singularities 1^ . Such 
pathologies could, in principle, prevent a numerical solver 
from evolving the spacetime in or near the critical regime. 

Now, as Garflnkle and Gundlach have discussed in de¬ 
tail [37|, an ideal coordinate system for numerical stud¬ 
ies of type II collapse is one which adapts itself to the 
self-similarity: for the DSS case this means that the 
metric coefficients and relevant matter variables are ex¬ 
actly periodic in the coordinates in the fashion given 
by (El). Glearly, if the coordinate system is adapted, then 
other than at the naked singularity—which is inacces¬ 
sible via finite-precision calculations—it should remain 
non-singular during a numerical evolution. One can then 
argue that ensuring that the numerical scheme has ade¬ 
quate resolution will be the key to successful simulation 
of the critical behaviour. At the same time, it is also clear 
that there will be coordinate systems which do not neces¬ 
sarily adapt but which nonetheless remain non-singular 
during critical collapse, at least over some range of scales, 
and which are therefore potentially useful for numerical 
calculations. We will see below that there is strong evi¬ 
dence that the coordinate systems we have used belong to 
the latter class, and weaker evidence that they do adapt 
to the self-similarity. 

Another potential source of problems, which is not spe¬ 
cific to hyperbolic formulations, relates to our restriction 
to spherical symmetry. As is well known, the singular 
points of curvilinear coordinate systems, r = 0 in our 
case, can sometimes require special treatment to ensure 
that numerical solutions remain regular there. In critical 
collapse the highly dynamical nature of the solution near 
r = 0 might naturally be expected to exacerbate prob¬ 
lems with regularity. In the work described below we 
have paid special attention to the ability of our approach 


to both fully resolve the near-critical configuration and 
maintain regularity of the solution at the origin. 

The remainder of this paper is organized as follows: in 
Sec.ini we review the generalized BSSN formulation and 
display the equations of motion for our model system. 
Sec. Uni expands the discussion of the issue of regularity 
at the coordinate singularity point, describes the numer¬ 
ical approach we have adopted, and provides details con¬ 
cerning the various tests and diagnostics we have used 
to validate our implementation. In Sec. m we present 
results computed using two distinct choices for the shift 
vector and provide conclusive evidence that the gener¬ 
alized BSSN formulation is capable of evolving in the 
critical regime in both cases. Sec. IVl contains some brief 
concluding remarks, and further details concerning the 
BSSN formalism in spherical symmetry and the scalar 
field equations of motion are included in App. and 
App. |B1 respectively. We adopt units where the gravi¬ 
tational constant and the speed of light are both unity: 
G = c = l. 


II. EQUATIONS OF MOTION 

The dynamical system we intend to study in the criti¬ 
cal collapse regime is a real, massless scalar field, self 
gravitating via Einstein’s equations, 

= 8ttT^, . (3) 

Here, is the energy-momentum tensor associated 
with the minimally coupled dt: 

Tf,, = , (4) 

and the evolution of the scalar field is given by 

V'"V^4' = 0. (5) 

The time-development of the geometry is then given by 
recasting Einstein’s equations as an evolution system 
based on the usual 3-1-1 expression for the spacetime met¬ 
ric: 

ds^ = —a^dt^ + "fij{dx^ + fi^dt){dx^ + f3^dt). ( 6 ) 

Here, the 3-metric components, 7 ij, are viewed as the 
fundamental dynamical geometrical variables and the 
lapse function, a, and shift vector, /3*, which encode the 
coordinate freedom of general relativity, must in general 
be prescribed independently of the equations of motion. 


A. Generalized BSSN 

We now summarize the BSSN formulation of Einstein’s 
equations and describe how it can be adapted to curvilin¬ 
ear coordinates. Readers interested in additional details 
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are directed to [s^ for a more pedagogical discussion. 

In the standard ADM formulation [3, the dy¬ 
namical Einstein equations are rewritten as evolution 
equations for the 3-metric and the extrinsic curvature 
{'yij,Kij}. The first difference between the BSSN for¬ 
mulation and the ADM decomposition is the conformal 
re-scaling of the ADM dynamical variables: 

lij = , (7) 

Kij = , ( 8 ) 

where e'^ is the conformal factor, 7 ^ is the conformal 
metric, Aij is the conformally rescaled trace-free part of 
the extrinsic curvature and K = Kij is the trace of 
the extrinsic curvature. Here by fixing the trace of Aij, 
and the determinant of the conformal metric, the set of 
primary ADM dynamical variables transforms to the new 
set: 

{77 ) lij ^ -^ij } ) ( 9 ) 

in the BSSN formulation. 


define 

Afc = yd (f) = f, ( 12 ) 

where denotes the Christoffel symbols associated with 
the flat metric. This definition makes this so-called con¬ 
formal connection. A*, a true vector and it becomes a 
primary dynamical variable in G-BSSN. 

We now summarize the G-BSSN equations, referring 
the reader to for more details, including a full deriva¬ 
tion. We begin by defining d ±, the time derivative oper¬ 
ator acting normally to the t = const, slices: 

d^ = dt-^0, (13) 

where denotes the Lie derivative along /3. We then 
have 

=-laK + , (14) 

6 6 

d±% = -2aAij - a‘^A,jDkl3^ , (15) 


In the original BSSN approach, the conformal met¬ 
ric is taken to have determinant 7 = 1 . However 
this choice is only suitable when we adopt coordinates in 
which the determinant of the flat-space metric reduces to 
unity. This is the case, of course, for Cartesian coordi¬ 
nates but is not so for general curvilinear systems. For 
instance, the flat 3-metric in spherical coordinates: 


d^K = --/^WjD,a+a{A,jA^^ + \K^)+4Tr{p+S ), (16) 
d±Aij = [—DiDja + a{Rij — SiTSij)]^^ 

+ a{KA,j-2AaA\)-a‘^A,jbkp'^, (17) 


ds^ = dr^ + r^de^ + sin^ OdijA , (10) 

has determinant 7 = sin^ 6. Recently, Brown [13 has 
resolved this issue by introducing a covariant version of 
the BSSN equations—the so-called generalized BSSN for¬ 
mulation, which we will hereafter refer to as G-BSSN—in 
which the primary dynamical variables are tensors so that 
the formulation can be adapted to non-Cartesian coordi¬ 
nate systems. In G-BSSN we no longer assume that the 
conformal 3-metric has determinant one. Rather, (j) be¬ 
comes a true scalar and for its dynamics to be determined 
a prescription for the time evolution of the determinant 
of jij must be given. In the following this will be done 
by requiring that the determinant be constant in time. 

Another main difference between the ADM decomposi¬ 
tion and BSSN is that the mixed spatial derivative terms 
occurring in the 3-Ricci tensor are eliminated through 
the definition of a new quantity, T^: 

ffc = , (11) 


= -2A^^d,a + j^^djdiP^ 


+ 2a 1 - -f^djK + eA^^djcj) 


-f 


3 L 




— 8717*-’ Sj 
(18) 


Here, a superscript TF denotes the trace-free part (with 
respect to the 3-metric 'jij) of a tensor, and Di is the co¬ 
variant derivative associated with the conformal metric 
jij. Additionally, the quantity a is an adjustable param¬ 
eter that is discussed below and typically is either 0 or 
1. Note that all the Lie derivatives in the G-BSSN equa¬ 
tions operate on true tensors and vectors of weight 0. For 
instance, 

^^A,j = . (19) 

Furthermore, in G-BSSN, rather than evolving dT51) . the 
redefined conformal connection. A®, is evolved via 


which becomes an additional, independent dynamical 
variable. Note that T* is not a vector as it is coordi¬ 
nate dependent. To extend this redefinition so that it 
is well suited for all coordinate choices, in G-BSSN we 


dtA’^ = , ( 20 ) 

where the time derivative 9 ( 7 *-^ is eliminated using (1151) . 
In equation CZl), Rij denotes the 3-Ricci tensor associ- 
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ated with 7 ^ and can be written as the sum 

Rij = Rf^ + Rij , ( 21 ) 


where Rfj is given by 

( 22 ) 

and Rij is the 3-Ricci tensor associated with the confor¬ 
mal metric: 

R^J = - if)f" + 


+ f-(2fff )fc^+fLffcz,) . (23) 

The matter fields p, S, Si and are defined by 

p = , (24) 

S = (25) 

, (26) 

S^J = , (27) 

where is the unit normal vector to the t = const, 
slices. 


As mentioned previously, we need to prescribe dynam¬ 
ics for the determinant of 7 ij to have a complete set of 
equations of motion for the G-BSSN dynamical variables. 
One approach is to hx the determinant to its initial value 
by demanding that 


da = Q. (28) 

This is the so-called Lagrangian option and is associated 
with the choice cr = 1 in the equations. Another option is 
to define the determinant to be constant along the normal 
direction to the time slices, which can be implemented 
by requiring c}j _7 = 0. This is usually referred to as 
the Lorentzian option, and is associated with the choice 
cr = 0. Here we choose (1^ . i. e. cr = 1. 

Note that in the G-BSSN equations the divergence of 
the shift vector, 


Dfc/3'= = ^5fc(v^^f, (29) 

no longer reduces to dkP^ since the determinant of the 
conformal metric ^ij is not necessarily 1 , but by virtue of 
the choice ( 1 ^ is equal to that of the initial background 
flat metric in the chosen curvilinear coordinates. 

As usual, when setting initial data for any given evo¬ 
lution of the coupled Einstein-matter equations we must 
solve the Hamiltonian and momentum constraints. In 


terms of the G-BSSN variables these are 

_ _ p<t> ^ _ 

^ = S^WiD.e* - R +— 

/ r J g g 

p 54 ’ 

- —K^ + 2ne^*p = 0, (30) 

= bj - ^e^^b^K - 87re®'^5* = 0 . (31) 


B. G-BSSN in spherical symmetry and gauge 
choices 


In spherical symmetry a generic form of the conformal 
metric jij is given by 


Jrrit,r) 0 0 

7ii = I 0 r^j00{t,r) 0 | . (32) 

0 0 r^ 7 ge(t, r) sin^ 0 

Similarly, a suitable ansatz for the traceless extrinsic cur¬ 
vature is 

/ Arr{t,r) 0 0 

Aij = I 0 r‘^Agg{t,r) 0 

\ 0 0 r‘^A00{t,r)sm^ 6 

(33) 

The shift vector and A* have only radial components: 

/3* = [/3(t,r),0,0] , (34) 


A* = 


A(t,r),0,0 


(35) 


Given (I32II35I) . the G-BSSN equations become a set of 
first order evolution equations for the 7 primary variables 

|())(t, r), 7rr-(i, r), ^gg{t, r),K{t, r), 
Arrit,r),Agg{t,r),A{t,r 


These are coupled to the evolution equation ([5]) for the 
scalar field and constrained by the initial conditions (1301 - 
I3I1). The explicit expressions for the full set of equations 
of motion are given in App. 

To fix the time slicing we implement a non-advective^ 


^ The terminology non-advective derives from the absence of 
an “advective” term, (d^dj, on the left hand side of equations 
(I36I3SD . we note that we also experimented with the advective 
versions of the equations, the results were very similar to those 
for the non-advective case; in particular, near-critical solutions 
exhibiting echoing and scaling could also be obtained. 
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version of the 1+log slicing condition:^ 


dta = —2ak . (36) 

for the spatial coordinates we either choose a zero shift: 

/3* = 0 , (37) 

or use what we will term the gamma-driver condition: 

dtP^ = /rA* - ry/3*. (38) 

Here, /r and r] are adjustable parameters which we set to 
= 3/4 and rj ~ 1/(2Madm), where Madm is the total 
mass of the system measured at infinity (see Sec. lIIIDTI) . 
We emphasize that ([38|) is not the usual Gamma-driver 
equation used in the standard BSSN approach: 

dtP^ = - r]p^, (39) 

but since it is a natural extension of the above to the G- 
BSSN case we have opted to use the same nomenclature. 
In the rest of this paper, we frequently refer to the shift 
vector evolved via (1551) as P'^. Explicitly, in spherical 
symmetry P’^ is defined by 

dtP^it,r) = nA{t,r) - r]P^(t,r). (40) 


III. NUMERICS 


We use a second order finite differencing method to 
discretize equations (I14II17I) and (l20l) . Further, the equa¬ 
tions of motion are transformed to a compactified radial 
coordinate that we denote by f and which is defined in 
terms of the original coordinate r by 




Rn 


r = e — e 


Rco 7 * Rqo ^ 


(41) 


where 5 and Roo are parameters with typical values 
i5 ~ —12 and Roo ~ 3. It is straightforward to verify 
the following: (1) the radial domain r = (0, oo) maps to 
the computational domain r = {S,Roo), (2) the deriva¬ 
tive dr/dr decreases toward the origin (f ~ 5), so that 
a uniform grid on f is a non-uniform grid on r with ap¬ 
proximately 10^ times more resolution close to the origin 
relative to the outer portion of the solution domain, f ~ 2 
(r ~ 10), where the support of the scalar held is initially 
concentrated, (3) the parameter 5 can be used to adjust 
the resolution near the origin; specihcally, decreasing 5 


increases the resolution near r = 0. For notational sim¬ 
plicity, however, in the following we omit the explicit 
dependence of the helds on r and denote the spacetime 
dependence of any dynamical variable X as previously: 
X (t,r(f)) = X{t,r). 

We use a hnite difference grid that is uniform in f 
and analytically transform all r-derivative terms in the 
equations of motion to their f-coordinate counterparts 
prior to hnite-differencing. 

We also developed a Maple-based toolkit [4l[ that au¬ 
tomates the process of discretizing an arbitrary derivative 
expression. This toolkit handles boundary conditions 
and generates a point-wise Newton-Gauss-Seidel solver 
in the form of Fortran routines for a given set of time de¬ 
pendent or elliptic PDFs . The calculations in this paper 
were all carried out using this infrastructure. 


A. Initialization 


The matter content is set by initializing the scalar held 
to a localized Gaussian shell: 

^(0,r) =pexp ^ , (42) 

where p, r^ and are parameters. Note that here r is 
the non-compactihed radial coordinate which is related to 
the compactihed coordinate f via dm. A typical initial 
prohle for the scalar held in our calculations has a,. ~ 
1, ro ~ 10, and p of order 10“^. We use the overall 
amplitude factor p as the tuning parameter to hnd critical 
solutions. We initialize the conformal metric ((32l) to the 
hat metric in spherical symmetry, 

>r(0,r) = 7ee(0,r) = 1, (43) 

and initialize the lapse function to unity, 

a(0,r) = l. (44) 

We also demand that the initial data be time-symmetric, 
A,(0,r) = iee(0,r) = A(0,r)=0, (45) 

/3(0,r) = A(0,r) = 0, (46) 


^ The reader can easily check that in the case of zero shift, 
the lapse choice given by ll36l l combined with m implies 
dt{(y- — ^/12) = 0. in Cartesian coordinates ^/12 = ln 7 , so 
this last equation gives a — In 7 = c(:r), where the function c(:r) 
is time independent, the choice c{x) = 1 then yields an algebraic 
expression for the lapse, a = 1 + ln 7 , which is the origin of the 
terminology “l+log slicing”. 


9t4'(t,r)|t=o = 0, (47) 

which means that the momentum constraint dm is 
trivially satisfied. This leaves the Hamiltonian con¬ 
straint (I30p which is solved as a two-point boundary value 
problem for the conformal factor at the initial time, 

V'(r) = . 


(48) 
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The outer boundary condition for ip, 

'0(r)|r=oo = 1 , (49) 

follows from asymptotic flatness, while at r = 0 we have 

drtpir)\r=o = 0 (50) 

since '0(*") must be an even function in r for regularity at 
the origin. 


B. Boundary Conditions 

Due to the fact that the metric has to be conformally 
flat at the origin we have 

Jrr{t,0) = J00{t,o). (51) 

Further, since we are using the Lagrangian choice, cr = 1, 
the determinant of jij must at all times be equal to its 
value at the initial time, so 

Irrlee = 1 • (52) 

From these two results we have 

>r-(t,0) = 7ee(i,0) = 1. (53) 

Using dsai) and (USD it is then easy to see that we must 
also have 

Arr{t,Q) = A0e{t,Q) = Q. (54) 

As is usual when working in spherical coordinates, 
many of the boundary conditions that must be applied at 
r = 0 follow from the demand that the solution be reg¬ 
ular there. Essentially, the various dynamical variables 
must have either even or odd “parity” with respect to 
expansion in r as r —>■ 0. Variables with even parity, typ¬ 
ically scalars or diagonal components of rank-2 tensors, 
must have vanishing radial derivative at r = 0, while odd 
parity functions, typically radial components of vectors, 
will themselves vanish at the origin. 

Applying these considerations to our set of unknowns 
we find 

drlrr{t,r)\r=o = dr^0e{t,r)\r=o = 0, (55) 

drArr{t,r)\r=0 = drA00{t,r)\r=O =0, (56) 

;9(t,0) = A(t,0)=0, (57) 

drK{t,r)\r=o = dr(t){t,r)\r =0 = l9r4'(t,»')|r=0 =0. (58) 

We use equations (|53I54I57D to fix the values of the func¬ 
tions at the origin and a forward finite-differencing of 


(l58l) to update K, 0 and 4' at r = 0. Further, we apply 
a forward finite-differencing of (155156^ to update the val¬ 
ues of the function at the grid point next to the origin. 
The H-log condition (I5S1) can be used directly at r = 0. 
Again, we emphasize that all of the r-derivative terms 
of the boundary conditions described above are analyti¬ 
cally transformed to the numerical coordinate, f, before 
the equations are finite-differenced. 

Since we are using compactified coordinates, all the 
variables are set to their flat spacetime values at the outer 
boundary r = oo: 

Irr =^00 = e'^ = a = \ at : (t, oo), (59) 

Ayr = A 00 = A' = A = /3 = d> = 0 at:(t, 00 ). (60) 

Here, we emphasize that spatial infinity, r = 00, corre¬ 
sponds to the finite compactified (computational) coor¬ 
dinate point f = i?oo- 


C. Evolution scheme and regularity 

We implemented a fully implicit, Crank-Nicolson 
finite differencing scheme to evolve the system of G-BSSN 
equations. The precise form of the continuum equations 
used is given in Add. [Aland all derivatives, both tempo¬ 
ral and spatial, were approximated using second-order- 
accurate finite difference expressions. 

During an evolution the correct limiting behaviour of 
the spatial metric components must be maintained near 
r = 0 to ensure a regular solution. For example, the 
limiting values of the conformal metric components ^rr 
and ^00 are given by 


^rr{t,r) = 1 -b O(r^), 

(61) 

I 00 {t,r) = 1-1- O(r^). 

(62) 


If the discrete approximations of the metric functions do 
not satisfy these conditions, then irregularity will mani¬ 
fest itself in the divergence of various expressions such as 
the Ricci tensor component (IA12I) 

i?rr = 2^^2_Z^ + ... , (63) 

r^lee 

which should converge to a finite value at the origin if 
conditions (|61I62D hold. 

One approach to resolve potential regularity issues is 
to regularize the equations [3, SO, Hi , by redefining the 
primary evolution variables, so that the equations be¬ 
come manifestly regular at the origin. Anothe r ap proach 
is to use implicit or partially implicit methods 1441. As re¬ 
cently shown by Montero and Cordeo-Carrion [45l| . such 
schemes can yield stable evolution without need for ex- 











plicit regularization. Baumgarte et al. also adopted a 
similar approach—using a partially implicit scheme with¬ 
out regularization—in an implementation of the G-BSSN 
formulation in spherical polar coordinates. 

As mentioned, our implementation is fully implicit and 
we have also found that our generalized BSSN equations 
can be evolved without any need for regularization at 
the origin, even in strong gravity scenarios where the 
spacetime metric has significant deviations from flatness 
near the origin. 

That said, we also experimented with other techniques 
aimed at improving regularity. For example, using the 
constraint equation dSH) and the fact that Aij is trace- 
free. 


^ -H 2^ = 0 , (64) 

Irr lee 

we can compute lee and Agg in terms of irr and Arr, 
respectively, rather than evolving them. However, when 
we did this we found no significant improvement in reg¬ 
ularity relative to the original scheme. 

Finally, to ensure our solutions remain smooth on the 
scale of the mesh we use fourth order Kreiss-Oliger dis¬ 
sipation in the numerical solution updates. 


D. Tests 

This section documents various tests we have made to 
validate the correctness of our numerical solver as well as 
the consistency of the finite-differencing method used to 
evolve the system of G-BSSN equations. We use a variety 
of diagnostic tools, including monitoring of the constraint 
equations, convergence tests of the primary dynamical 
variables, and a direct computation to check if the metric 
and matter fields calculated via the G-BSSN formulation 
satisfy the covariant form of Einstein’s equations. All of 
the calculations were performed using the 1-1-log slicing 
condition ((Ml) and either /3 = 0 or /3 = where /3® 
satisfies the Gamma-driver condition (HUl) . 


1. Constraints and conserved quantities 

We monitor the evolution of the constraint equations 
(I30I31|) during a strongly-gravitating evolution where the 
nonlinearities of the equations are most pronounced. As 
demonstrated in Fig.[T](a), at resolutions typical of those 
used in our study, the Hamiltonian constraint is well pre¬ 
served during such an evolution and, importantly, the 
deviations from conservation converge to zero at second 
order in the mesh spacing as expected. 

The total mass-content of the spacetime seen at spatial 
infinity (the ADM mass) is a conserved quantity. Here, 
using the G-BSSN variables the Misner-Sharp mass func- 



t/M 





t/M 


FIG. 1: Results from various tests that verify the accuracy 
and consistency of our numerical solver and the finite dif¬ 
ferencing method used to integrate the equations, (a) The 
evolution of the / 2 -norm (RMS value) of the Hamiltonian con¬ 
straint. The norm is plotted for 3 different resolutions h, h/2 
and /i/4 corresponding to Nr = 512, 1024 and 2048, respec¬ 
tively. The data for the Nr = 1024 and Nr = 2048 curves 
have been rescaled by factors of 4 and 16, respectively, and 
the overlap of the three lines thus signals the expected sec¬ 
ond order convergence to zero of the constraint deviation. 
We observe similar convergence properties for the momentum 
constraint as well as the constraint equation (I12II for A*, and 
the constraint (1641) for the trace of Aij. Additionally, since 
the operator used to evaluate the residual of the Hamiltonian 
constraint is distinct from that used in the determination of 
the initial data, the test also validates the initial data solver, 
(b) Conservation of the ADM mass during the evolution of 
strong initial data. Here the deviation of the mass from its 
time average is plotted for 3 different resolutions. Higher res¬ 
olution values have again been rescaled so that overlap of the 
curves demonstrates 0{h^) convergence to 0 of the deviation 
of the total mass, (c) The convergence factor defined in (1681) 
for three of the primary BSSN variables: 7rr, K, and Age- In 
the limit h —>■ 0 we expect all curves to tend to the constant 
4. The plot thus provides evidence for second order conver¬ 
gence of all of the values throughout the evolution. All of the 
other primary BSSN variables as well as the dynamical scalar 
field quantities demonstrate the same convergence, (d) Direct 
verification that the metric found by numerically solving the 
BSSN equations satisfies Einstein’s equations in their covari¬ 
ant form. Here the tr component of the residual defined in 
([72)) is plotted for 3 different resolutions. Once more, higher 
resolution values have been rescaled so that overlap of the 
curves signals the expected 0{h^) convergence of the resid¬ 
uals to 0. All of the plots correspond to evolution of strong 
subcritical initial data with H-log slicing. For (a) and (b) the 
shift vector was set to 0, while in (c) and (d) it was evolved 
using the Gamma-driver condition (i.e. /3 = 
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tion is given by 


Mir) = 





2 


^ 2-f0e e'^ J 


(65) 

The total mass, Madm, can be evaluated at the outer 
boundary, 


-Madm = M(r = oo). (66) 

The deviation of the total mass from its time average is 
plotted in Fig. HDb); as the resolution of the numerical 
grid increases the variations converge to zero in a second 
order fashion. 


2. Convergence test 

As mentioned in Sec. IIII Al and Sec. IIIICl we imple¬ 
mented our code using second order finite differencing of 
all spatial and temporal derivatives. Denoting any con¬ 
tinuum solution component by q(t,X), where X is the 
spatial coordinate, and a discrete approximation to it 
computed at finite difference resolution, h, by q^{t,X), 
to leading order in h we expect 

q'^it,X) = q{t,X) + h‘^e 2 [q]it,X) + ... . (67) 

Fixing initial data, we perform a sequence of calculations 
with resolutions h, h/2 and h/A and then compute a con¬ 
vergence factor, C{t;q), defined by 


using the primary G-BSSN variables, 7y and (j). In par¬ 
ticular, a and b are simply given by 


■iS 

II 

(70) 

b{t,r) = e^‘^‘^*’’'^ee(f,r). 

(71) 


We then check to see if the metric (IMl) satisfies the co¬ 
variant Einstein equations m to the expected level of 
truncation error. Specifically, defining the residual 


= G't - , (72) 


and replacing all derivatives in with second order fi¬ 
nite differences, we expect to converge to zero as 
0(h?) as /i —?> 0.“^ Precisely this behaviour is shown in 
Fig. Ed). This is a particularly robust test of our imple¬ 
mentation since the non-trivial components of the covari¬ 
ant Einstein equations are quite complicated and, super¬ 
ficially at least, algebraically independent of the BSSN 
equations. For instance, the tr component of the resid¬ 
ual (1721) is given by 


E\ 


2/3 /dra ^drb 9^0\ 
ra^ \ a b a J 
2/3 / d^b drttdrb {dr'A>)'^ dradrb\ 

\ b ab 2 ab ) 


2 /dtdrb drAfdtA! dtadrb dtbdra\ 
\ 6 2 ab ab J 


ra^ \ a b J 


(73) 


\\q\t,X)-q^/\t,X )\\2 


( 68 ) 


where 11-112 is the I 2 norm, i.e. the root mean square 
(RMS) value. It is straightforward to argue from (l67l) 
that, for sufficiently small /i, C{t,q) should approach 4 if 
the solution is converging at second order. The values of 
the convergence factor for a selection of dynamical vari¬ 
ables are plotted for a strong-data evolution in Fig. [TKc) 
and provide clear evidence that the solution is second- 
order convergent throughout the time evolution. 


3. Direct validation via Einstein’s equations 


A direct method to test the fidelity of our numerical 
solver involves the evaluation of a residual based on the 
covariant form of Einstein’s equations. We start with a 
reconstruction of the four-dimensional metric in spherical 
symmetry, 

= (—-I- /3^a^)dt^ -I- 2a^f3dtdr + -|- r'^b'^dG? , 

(69) 


and depends on all of the dynamical variables of the sys¬ 
tem. The observed convergence of the residual is only 
plausible if (I) our G-BSSN equations (11411181) have been 
correctly derived from the covariant Einstein equations, 
(2) we have discretized the geometric and matter equa¬ 
tions properly, and (3) we have solved the full set of dis¬ 
cretized equations correctly. 


E. Finding black hole threshold solutions 

The strength of the initial data can be set by adjust¬ 
ing the amplitude of the scalar field, p, in (H^ . For weak 
enough initial data (small enough p), the matter shell will 
reach the origin and then disperse, with the final state 
being a flat spacetime geometry. Sufficiently strong ini¬ 
tial data, on the other hand (large enough p), results in a 
matter concentration in the vicinity of the origin which is 
sufficiently self-gravitating that a black hole forms. Us- 


Although it is not crucial for the usefulness of this test, we dis¬ 
cretize the Eii using a difference scheme that is distinct from the 
one used in the main code. 
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ing a binary search, we can find the threshold initial data, 
defined hy p = p* , for which p < p* results in dispersal 
while p > p* yields black hole formation. At any stage 
of the calculation, the binary search is defined by two 
“bracketing” values, pi and ph, such that evolutions with 
p = pi and p = Ph result in dispersal and black hole 
formation, respectively. It is convenient to define the 
amount of parameter tuning that has occurred by the 
dimensionless quantity 

5p = Ph-Pi ( 74 ) 

Pi 


The dispersal case can be detected easily as the scalar 
field leaves the vicinity of the origin and the geometry 
approaches flat spacetime. To detect black hole forma¬ 
tion, we use an apparent horizon finder to locate a surface 
r = const, on which the divergence of the outgoing null 
rays vanishes. We first define the divergence function 

0 = , (75) 


where is the induced metric on the constant r surface. 
In spherical symmetry with metric (I69p we have 

= diag (0,0,r^6^,r^6^sin^6») , (76) 

where is the null outgoing vector given by 


kfj, = —= [aj3 - a, a, 0, 0] . 

v2 


(77) 


Therefore, dZB becomes 

^ = ^{-dt{h)+{--^]dr{rb) 
ro \a \a a 


(78) 


The formation of an apparent horizon^ is signaled by 
the value of the function 0 crossing zero at some radius 
and, modulo cosmic censorship, implies that the space- 
time contains a black hole. We note that since the focus 
of our work was on the critical (threshold) solution we 
made no effort to continue evolutions beyond the detec¬ 
tion of trapped surfaces. 



ln(T,-T) 


FIG. 2: Echoing behaviour in the scalar field for a marginally 
subcritical evolution with 5p « 10“*^^. The main plot displays 
the central value of the scalar field versus a logarithmically 
scaled time parameter, ln(r/ — T), where T is central proper 
time and T/ is the approximate value of that time when near- 
critical evolution ceases and the total dispersal of the pulse 
to infinity begins. This particular scaling is chosen solely to 
more clearly demonstrate the evolution of the central value of 
during the critical phase through to dispersal. Note that 
our choice of abscissa means that evolution proceeds from 
right to left. The inset also plots at r = 0 but now in the 
“natural” logarithmic time coordinate r = ln(r* — T) where 
T* is the “accumulation time” at which the solution becomes 
singular and which has been estimated based on the posi¬ 
tions of the extrema in The amplitude of the scalar field 
at the origin oscillates between (—0.61,0.61), consistent with 
the calculations reported in [T^. The data yield an echoing 
exponent of A = 3.43 ± 0.02 which is in agreement with the 
value A = 3.445452402(3) Martin-Garcia and Gundlach have 
computed by treating the computation of the precisely-critical 
solution as an eigenvalue problem (^. 


A. Zero shift 


IV. RESULTS 

In this section we describe results from two sets of nu¬ 
merical experiments to study the efficacy of the G-BSSN 
formulation in the context of critical collapse. Again, our 
calculations use the standard I-blog slicing condition for 
the lapse, and a shift which is either zero or determined 
from the Gamma-driver condition. 


® Technically a marginally trapped surface—the apparent horizon 
being the outermost of these. 


We first perform a collection of numerical experiments 
where the shift vector is set to zero. As described in 
Sec. IIII El in principle we can find the black hole thresh¬ 
old solution p ~ p* using a binary search algorithm which 
at any stage is defined by two values pi and ph, with 
Pi < p < Ph, and where pi corresponds to dispersal 
(weak data) while ph corresponds to black hole forma¬ 
tion (strong data). 

As discussed in the introduction, the massless scalar 
collapse model has a very well-known critical solution, 
and we summarize the features most relevant to our 
study here. The threshold configuration is discretely self¬ 
similar with an echoing exponent measured from the first 
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ln|p-p‘| 

FIG. 3: The msiximum central value, TZmax, of the four¬ 
dimensional Ricci scalar, TZ, attained during subcritical evo¬ 
lution as a function of the logarithmic distance In \p — p*\ of 
the tuning parameter from the critical value. As first ob¬ 
served in the Ricci scalar scales as i? ~ \p — p*\~^'^, where 
7 is the universal mass-scaling exponent in ©. The value 
7 = 0.38 ± 0.01 computed via a least squares ht is in agree¬ 
ment with the original calculations [l^ as well as many other 
subsequent computations. We note that the oscillations of 
the data about the linear ht are almost certainly genuine, at 
least in part. As discussed in the text, we expect a periodic 
wiggle in the data with period A/( 27 ) ~ 4.61. Performing 
a Fourier analysis of the residuals to the linear ht we hnd a 
peak at about 4 with a bandwidth of approximately 1, con¬ 
sistent with that expectation. As described in more detail 
in the text, although we have data from computations with 
In Ip — p*| < —25, we do not include it in the ht. The naive 
method we use to estimate p* means that the relative uncer¬ 
tainty in p—p* grows substantially as p —^ p* so that inclusion 
of the data from the most nearly-critical calculations will cor¬ 
rupt the overall ht. 


calculations to be A Ri 3.44 [T^. Following the origi¬ 
nal studies, Gundlach [l^ showed that the construction 
of the precisely discretely self-similar spacetime could be 
posed as an eigenvalue problem, the solution of which led 
to the more accurate value A = 3.4439 ±0.0004. This es¬ 
timate was subsequently improved by Martin-Garcia and 
Gundlach to A = 3.445452402(3) [H. 

The original calculations determined a value 7 r; 0.37 
for the mass-scaling exponent ; further work based on 
perturbation theory gave 7 r: 0.374 [H, Here it is 
important to note that, as pointed out independently by 
Gundlach [l^ and Hod and Piran the simple power 
law scaling o gets modified for discretely self-similar 
critical solutions to 

In M = 7 In Ip — p* I ± c ± / (7 In Ip — p* I ± c) , (79) 



r 


FIG. 4: Discrete self-similarity of the geometry of spacetime 
in the black hole threshold evolution previously discussed in 
Fig. O Here, the G-BSSN variable cj) is plotted as a function 
of the computational radial coordinate f at the accumulation 
time t*. Note that from (pl) measures the deviation of 
the determinant of the 3-metric from that of a flat metric. 
The inset graph is the radial derivative of (f) scaled by y/r to 
highlight the formation of fine structure in the geometry of 
the critical solution. The approximate periodicity of y/r4>' in 
ln(r) (modulo an overall varying scale) provides weak evidence 
that the coordinate system used in the calculation adapts to 
the self-similarity of the critical solution. 


where / is a universal function with period A and c is 
a constant depending on the initial data. This results in 
the superposition of a periodic “wiggle” in the otherwise 
linear scaling of In M as a function of In |p — p* |. 

Finally, Garfinkle and Duncan [2lj pointed out that 
near-critical scaling is seen in physical quantities other 
than the mass and, dependent on the quantity, in the 
subcritical as well as supercritical regime. In particular 
they argued that in subcritical evolutions the maximum 
central value, T^max, of the four-dimensional Ricci scalar, 
TZ, defined by 

7^^lax = max7^(0, t ), (80) 

should satisfy the scaling 

(81) 

where the factor —2 in the scaling exponent can be 
deduced from the fact that the curvature has units of 
length”^. For the discretely self-similar case this scaling 
law is also modulated by a wiggle with period A/( 27 ), 
which for the massless scalar field is about 4.61. 

Using initial data given by (P|) we tune p so that it 
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FIG. 5: Snapshots of radial mass density for a marginally 
subcritical calculation {5p « 10“^^, Nr = 2048). Plotted is 
dm/dr = r^p{t,r) where p is defined by (I24II . In this calcu¬ 
lation = 0 so we also have dm/dr = T\. As the solution 
evolves, development of echos is clearly seen. In the hnal 
frame, which is at an instant t = 15.8 that is close to the 
accumulation time t*, we observe 4 echos. Note that we do 
not count the tall thin peak at the extreme left nor the hrst 
two peaks on the right as echos. The skinny peak will develop 
into an echo as p is tuned closer to p*. The two peaks on the 
right account for the bulk of the matter and represent the 
part of the initial pulse that implodes through the origin and 
then disperses “promptly”, i.e. without participating in the 
strongly self-gravitating dynamics. A corresponding plot for 
an evolution far from criticality would contain only those two 
peaks. Note that the hrst three plots use the computational 
coordinate f to provide a sense of the actual numerical calcu¬ 
lation, while the last plot uses ln(r) in order to best highlight 
the discrete self-similarity of the threshold solution. As is the 
case for the data plotted in the inset of the previous hgure, 
the approximate periodicity of the mass density in ln(r) sug¬ 
gests that the coordinates are adapting to the self-symmetry 
of the critical spacetime. 


is close to the critical value: typically this involves re¬ 
ducing the value of 5p defined by dZil) so that it is about 
10“^^, which is a few orders of magnitude larger than ma¬ 
chine precision. Our implementation includes code that 
actively monitors the dynamical variables for any indi¬ 
cations of coordinate singularities or other pathologies 
which could cause the numerical solver to fail. Provided 
that such pathologies do not develop, we expect to ob¬ 
serve features characteristic of critical collapse—discrete 
self-similarity and mass scaling in particular—to emerge 
as p ^ p*. 

One way the discrete self-similarity of the critical 


solution is manifested is as a sequence of “echoes” — 
oscillations of the scalar field near the origin such that 
after each oscillation the profile of the scalar field is re¬ 
peated but on a scale exp(A) smaller than that of the 
preceding echo (see Eq. ([2])). The oscillations are simi¬ 
larly periodic in the logarithmic time scale ln(T* — T), 
where T is the proper time measured at the origin, 

T{t)= [ a{i,0)di, (82) 

Jo 

and T* is the accumulation time at which the singular¬ 
ity forms (always at r = 0). Furthermore, viewed at the 
origin, the oscillations of the scalar field occur at a fixed 
amplitude of about 0.61 (with our units and conventions 
for the Einstein’s equations). As shown in Fig. [2J when 
we tune the initial data to the critical value, the central 
value of scalar field exhibits oscillatory behaviour and 
the amplitude is close to the expected value. The an¬ 
ticipated periodicity in logarithmic time is also apparent 
with a measured A = 3.43 ± 0.02, in agreement with 
previous results. We thus have strong evidence that the 
evolution has indeed approached the critical regime and 
that the measured oscillations are true echos rather than 
numerical artifacts. 

Evidence that our code correctly captures the expected 
critical scaling behaviour (ISTl) of TZmax is presented in 
Fig. [3] We find 7 = 0.38 ±0.01, consistent with pre¬ 
vious calculations. We note that we can measure scal¬ 
ing from our computations up to \n\p — p*\ = —29 (or 
\p—p*\ ~ 10“^^). However, in Fig. Owe have excluded the 
last few values closest to the critical point from both the 
plot and the linear fit: specifically, we truncate the fit at 
\a\p—p*\ = —25. The rationale for this is that we use the 
largest subcritical value of p as an approximation to the 
critical value p* rather than, for example, implementing a 
multi-parameter fit that includes p* as one of the param¬ 
eters. Our estimate of p* thus has an intrinsic error of 
PS 10“^^ and by fitting to data with In \p—p*\ > —25 
we render the error in the p* estimate essentially irrel¬ 
evant. We note that consistent with the early observa¬ 
tions of the robustness of mass scaling in the model [l^ , 
measuring the exponent 7 can be achieved by moder¬ 
ate tuning, in this case ln|p —p*| Ri —9, (i.e. 5p k, 10“^). 
However, to be able to observe the echoing exponent (the 
oscillations around the fitting line, for example) we need 
to tune much closer to the critical value. 

The echoing behaviour of the critical solution is also 
reflected in the geometry of spacetime and the matter 
variables other than the scalar field. Fig. 2] shows the 
radial profile of the G-BSSN variable </> at an instant 
close to the accumulation time T*. As seen in this plot, 
fine structure develops in the function in the near-critical 
regime. Observe that from the definition © and the 
choice (1^51) . the scalar (p is the ratio of the determinant 
of the 3-metric, 7, to the determinant of the flat metric. 
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(t>= (83) 

The radial matter density, dm/dr = r'^p, is a conve¬ 
nient diagnostic quantity for viewing near-critical evolu¬ 
tion. Snapshots of this fnnction from a typical marginally 
subcritical calculation are shown in Fig. [5j the echoing 
behavior is clearly evident in the sequence. The number 
of echos is dependent on the degree to which the solu¬ 
tion has been tuned to criticality. In this case, where 
Sp = 10“^^, we expect and see about 4 echos (last frame 
of the figure). Here we note that each of the echos in 
dm/dr corresponds to half of one of the scalar field oscil¬ 
lations shown in Fig. [5] (where the inset shows about 2^ 
full cycles). 

Fig. EK a) plots the central matter density p{t, 0) for 
a marginally supercritical calculation. In accord with 
the self-similar nature of the near-critical solution, the 
central density grows exponentially with time. Fig. EKb) 
is a snapshot of the extrinsic curvature at the critical time 
t K, t* while Fig. [61(c) shows the dynamics of the central 
value of the lapse function and compares it with a from 
the calculations performed with fd = jdc described in the 
next section. Fig.EKd) displays the profile of the lapse at 
the critical time. 

We note that we have not fully resolved the critical be¬ 
haviour in the sense of having tuned p to the limit of ma¬ 
chine precision, Sp ~ 10“^®, which would capture roughly 
2 additional echos in the threshold solution (one full echo 
in the scalar field). In principle, by setting Nr sufficiently 
large we could almost certainly do so since there are no 
indications that our method would break down at higher 
resolution and closer to criticality. However, we estimate 
that the required compute time for a complete critical 
search would increase from weeks to several months and 
we have thus not done this. Ultimately, a more effective 
approach to enhancing the resolution would be to incor¬ 
porate a technique such as adaptive mesh refinement into 
our solver. 

The results displayed in Figs. 2-6 provide strong ev¬ 
idence that the coordinate system consisting of l-f-log 
lapse and zero shift remains non-singular in the criti¬ 
cal regime, at least for the range of scales probed for 
Sp ~ 10“^^. Additionally, the approximate periodicity 
in ln(r) that can be seen, for example, in (Fig. 4) 

and dm/dr (Fig. 5) suggests that the coordinates may 
be adapting to the self-similarity. Whether or not this is 
actually the case is a matter requiring further study. 



FIG. 6: Profiles of matter and geometry variables from 

strongly gravitating, near-critical evolutions where the echo¬ 
ing behaviour emerges. Results were computed using U-log 
slicing and zero shift, except for the dashed curve in (c) where 
d = (a) Central energy density, p(T, 0), as a function 

of proper central time, T, and in logarithmic scale for a su¬ 
percritical evolution. The density oscillates and grows expo¬ 
nentially as the system approaches the critical solution and 
then eventually collapses to form a black hole, (b) Profile of 
the extrinsic curvature, K(t*,r )—scaled by in order to 
make the echoing behaviour more visible—where t* denotes 
a time very close to the accumulation time. The evolution is 
marginally subcritical in this case, (c) Central value of the 
lapse function, a, during subcritical evolutions with (5 = 0 
(solid) and d = do (dashed). The plots use a logarithmically 
transformed proper time variable, — ln(r/ — T), where Tf is 
the approximate time at which the final dispersal of the pulse 
from the origin begins. In both cases a exhibits echoing and 
there is no evidence of pathological behaviour, such as the 
lapse collapsing or becoming negative. The close agreement 
of a for the two choices of (5 indicates that the time slicing 
varies little between the two coordinate systems. Note that 
there are three extra oscillations for the /3 = Pc case, in the 
time interval — ln(r/ — T) > 4.5. These are spurious and 
due to a lack of finite-difference resolution; there are only 6 
time steps in each oscillation, (d) Radial profile a{t*,r) at a 
time t = t* which is close to the accumulation time and when 
the self-similarity and echoing in the spacetime geometry is 
apparent. 


B. Gamma-driver shift 

We now briefly report on experiments similar to those 
of the previous section but where the shift was evolved 
with the Gamma-driver condition psp . A principal ob¬ 
servation is that this gauge also facilitates near-critical 


evolutions with results similar to the /3 = 0 choice. In 
particular, we are again able to observe all of the char¬ 
acteristics of the black hole threshold solution. 

The gauge condition (1551) acts as a damping factor for 
the conformal connection. A*, and we would therefore ex¬ 
pect to observe a signihcant change in the profile of A* at 
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FIG. 7: Profiles of various G-BSSN variables from marginally 
subcritical evolutions, (a) Profile of the conformal connection 
A as computed with P = Q and at a time t* close to the ac¬ 
cumulation time. Note that the function has been scaled by 
r and in fact diverges like 1/r. (b) Profile of A as computed 
with P = P^, again at a moment close to the accumulation 
time. Here the 1/r growth seen when the condition /? = 0 
is adopted is absent, (c) Profile of the central spatial deriva¬ 
tive of the shift vector, P'{t,0), as computed with P = p'^. 
As the echos develop closer to the origin, P' increases and 
presumably will diverge in the continuum, precisely-critical 
limit, (d) Time development of the Z 2 -norm of the extrinsic 
curvature during subcritical evolutions for both the P = 0 
and P = P^ calculations. In both cases the extrinsic curva¬ 
ture develops a divergent profile near r = 0 in the critical 
regime. 


threshold relative to the zero-shift case. This expectation 
is borne out by the comparison illustrated in Figs. 0 (a) 
and (b). When /3 = 0, A* diverges as 1/r close to the ori¬ 
gin while it appears to have finite amplitude for /3 = . 

We find that the shift develops very sharp oscillations 
near the origin; some typical behaviour can be seen in 
the plot of P'{t,0) shown in Fig. [7ljc). We believe that 
these oscillations are genuine and our expectation is that 
P'(t,0) will diverge in the precise critical limit. Further, 
we observe that the oscillations can create numerical ar¬ 
tifacts and generally require higher resolution relative to 
the P = 0 case, as well as dissipation, to be controlled. 
Indeed, when using the Gamma-driver condition we find 
that Kreiss/Oliger dissipation is crucial to suppress un¬ 
resolved high frequency oscillations close to the origin. 
Fig. [3 (d) shows the growth in the norm of the extrinsic 
curvature during a subcritical evolution. The norm of 
K does not exhibit any significant difference for the two 


choices of the shift. 

As was the case for the P = 0 calculations, the results 
shown in Figs. 6 and 7 strongly suggest that the combi¬ 
nation of 1+log slicing and Gamma driver shift provides 
a coordinate system which is adequate for computing the 
near-critical solution. In addition, the approximate peri¬ 
odicity seen in Figs. 6(b), 6(c), 7(a) and 7(b) suggest that 
this gauge may also be adapting to the self-symmetry. 


V. CONCLUSION 

We have described a numerical code that implements a 
generalized BSSN formulation adapted to spherical sym¬ 
metry. Using standard dynamical coordinate choices, in¬ 
cluding 14-log slicing and a shift which either vanished 
or satisfied a Gamma-driver condition, we focused specif¬ 
ically on the applicability of the formulation and the 
gauge choices to studies of type II critical phenomena. 
As a test of the approach we revisited the model of mass¬ 
less scalar collapse, where the properties of the critical 
solution are very well known from previous work. For 
both choices of the shift, we found that our code was 
able to generate evolutions that were very close to crit¬ 
icality so that, in particular, we could observe the ex¬ 
pected discrete self-similarity of the critical solution. To 
our knowledge, this is the first fully evolutionary imple¬ 
mentation of a hyperbolic formulation of Einstein’s equa¬ 
tions that has been able to unequivocally resolve discrete 
self-similarity in type II collapse. Furthermore, mea¬ 
sured properties from near-critical solutions, including 
the mass-scaling and echoing exponents, are in agreement 
with previous work. Our results strongly suggest that 
the G-BSSN formulation, in conjunction with standard 
dynamical coordinate conditions, is capable of evolving 
the spacetime near criticality without the development of 
coordinate pathologies. There is also some evidence that 
both gauges adapt to the self-similarity, but we have not 
yet studied this issue in any detail. 

We found that certain of the primary G-BSSN vari¬ 
ables diverge as the critical solution is approached: this 
is only to be expected since the precisely critical solution 
contains a naked singularity. Dealing with such solu¬ 
tion features in a stable and accurate manner presents 
a challenge for any code and in our case we found that 
a combination of a non-uniform grid and Kreiss/Oliger 
dissipation was crucial. Our use of a time-implicit evolu¬ 
tion scheme may have also been important although we 
did not experiment with that aspect of our implementa¬ 
tion. However, we suspect that the implicit time-stepping 
helped maintain regularity of the solutions near r = 0, 
as other researchers have found. 

Given the success of the G-BSSN approach, it is natu¬ 
ral to consider its generalization and application to set¬ 
tings with less symmetry, but where curvilinear coordi¬ 
nates are still adopted. In particular, one axisymmetric 
problem that has yet to be resolved is the collapse of 
pure gravity waves. This scenario arguably provides the 
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most fundamental critical phenomena in gravity as the 
behaviour must be intrinsic to the Einstein equations, 
rather than being dependent on some matter source. 
Critical collapse of gravitational waves-with mass scal¬ 
ing and echoing—was observed by Abrahams and Evans 
over 20 years ago [l^. However, their original results 
have proven very difficult to reproduce (or refute) [28l - 
[m, HJ. We refer the reader to the recent paper by 
Hilditch et al. [s^ for detailed discussions concerning 
some apparent inconsistencies among the follow-up stud¬ 
ies, as well as the challenges and complications involved 
in evolving various types of nonlinear gravitational waves. 
We are currently extending the methodology described 
above to the axisymmetric case with plans to use the re¬ 
sulting code to study vacuum collapse. Results from this 
undertaking will be reported in a future paper. 


^ r. 196 r 

^99 = rdra- -h —Ora 

'Trr ^ 


f drl99 
\ 'Trr 



The trace of S^ij is 


^ = e-40 + 2 ^ 

V Irr r‘‘l99 

Then the evolution of K is given by 


dtK = + + 


^ Irr 

+ pdrK + 4:Tra{p + S) 


%9 


(A7) 


(A8) 


(A9) 


and the evolution equations for the traceless part of the 
extrinsic curvature are 
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Appendix A: BSSN in Spherical Symmetry 

In this appendix, we provide the explicit expressions of 
the G-BSSN evolution equations in spherical symmetry. 

The evolution equations (11411151) for (j) and the compo¬ 
nents of the conformal metric lij simplify to 


dt4> = \aK + ISdrcj) + crl‘ 
6 6 


(Al) 


dtlrr = —2aArr + jSdrlrr + 2%rdrl3 — CT-lrr^ , (A2) 

O 


(3 2 

dtlee = - 2 aAg 0 + pdriee + ‘ 2^—196 - , (A3) 

r 3 

where 33S is the divergence of the shift vector, 


^(t, r) = A/3* = ^ 

r V 2 'yrr lee 


(A4) 


To display the equation of motion for the trace of the 
extrinsic curvature K and Ay we hrst define 

= DiDja , (A5) 

which has 2 independent components, 

^rr = d^a — dra (4,dr(j^ , (A6) 

V Irr J 


dtArr = e + a [Rjf + SttSJ/')] 

+ a( ArrK - 


-I- 2Arrdrl3 + POrArr — a-£§Arr , (AlO) 

o 


+ a{Rjf -kSTrS'J/’)] 

+ a (AggK — 2^::^^ 

\ lee J 

+ 2 — Agg-\-pdrAgg — a — A 00 ^, (All) 

r 3 

where R denotes the 3-Ricci tensor with non-vanishing 
components 




i{drlrrf {drleef , ~ a a , ^ O ~ A 

- — - ^2 - -|- —Ur'^rr-^ 




‘^'^99 


^ (A ^ J. ^rlrr 2drl00 2^rrdrlee 

V lee T 00 


— 4:d^(j) + 2 dr 4 > 


drlrr drlBB 


Irr 199 
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(drP^ - 2 dlP - A{drPf 
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drigg 6dr4>lee 


lee 


Ir 
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(A13) 


In the above expressions the superscript TF denotes ap- 
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plication of the trace-free-part operator, whose action can 
be written explicitly as 

Xjf = Xrr - \lrrX = U Xrr - , (A14) 

3 3 V leer-^ ) 

Xjf = Xee - \le 0 X = i (xgs - . (A15) 

0 \ 'Trr / 


definition 

Sn = — - ^dr'I'R) . (B4) 

a 

We then find the following evolution equations for 
and Er: 

dt'^R = j^Sr +/Sdr'i'R , (B5) 


Here X represents any of the tensors 3), R 01 S. 
Finally, the evolution of A® reduces to 


9tA = fidrR — drfiR + t— 

'Jrr 

a I drjrr^rr 29 , 


^rr \ ^rr 3 , 


7r 


^2 

Irr 


- -+ 

nis 


,' 2 ~ dr^\ 2 

+ cr ( + —— ) H—z— 

. 3 3jrr J r-fgg 


9r/3 - 
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, drUAr 


^2 

Irr 
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dt^R 


ab"^ 

a 



^r + 2 


dr^R ^ 


+ dr^ rOt 


Pdr^R + ^RdrP + ^R - 

r 

aab^d'^^^Vm). 



(B6) 


The evolution equations for 4// and Ej follow from the 
index substitutions R ^ I in the right hand sides of (IB5I1 
and (IB6I) . respectively. 

The matter source terms in the G-BSSN equations, 
namely p, S, S'*, Sy, can be simplified by defining H and 
$ as 


Appendix B: Scalar field dynamics and 
energy-momentnm tensor in spherical symmetry 


Here we present the evolution equations of a complex 
scalar field, with an arbitrary potential V, minimally cou¬ 
pled to gravity. The governing equations for a massless 
real scalar field follow as a special case where the imag¬ 
inary part of the field and the potential are both set to 
zero. 

The geometry of spacetime is given by a generic metric 
in spherical symmetry: 

ds^ = (-a^ + /3^a^)dt^ + 2a^pdtdr + a^dr^ + r'^b^dn^, 

(Bl) 

where a, &, a and /3 are all functions of t and r and 
where a and b are related to the primary BSSN variables 
via a = ^rr exp(40) and b = ^gg exp(4(/)). 

The complex scalar field is given in terms of real and 
imaginary parts, 4'^ and respectively, 

^ ^'^Rit,r)+i^i{t,r), (B2) 


H = - (StT - Pdr'if) = UR{t, r) + ini{t, r ), (B7) 

a 


Ar = - {dt'fR - Pdr^R) = ^ , (B8) 

a 

B/ = - - Pdr-^i) = ^ , (B9) 

a 0 ^ 

^ = dr'i'= ^R{t,r) + , (BIO) 

<i>R = dr'fR, (Bll) 

$^ = 9^4'/. (B12) 


Using these definitions, the variables p and S are given 
by 


P{ix) 


|np + |4>p U(|4/|) 

2 a 2 2 


(B13) 


and has an associated energy-momentum tensor 

-f V^V^4'j - 

- \9>.uVm)- (B3) 


3|nP-|4>P 3 

S{t,r) = ' ' - -U(|4^|). (B14) 

In spherical symmetry, S'* has only a radial component, 
S* = [S’'(t,r),0,0] , (B15) 


with 


The evolution of the real part of the scalar field can be 
reduced to a pair of first-order-in-time equations via the 


gr 


-I- Hjtl’/ 
a 


(B16) 
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Similarly, the spatial stress tensor, Sij, has only two in¬ 
dependent components. 


Srr(t,r) 0 0 

s,, = ( 0 r^See 0 

0 0 sin^ 9S00 


(B17) 


with 


,B18, 


S00 = 


V 2a2 


vm) ^ 


(B19) 
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